INTRODUCTION 


Monte Carlo transport codes like FLUKA are useful for many purposes, and one 
of those is the simulation of the effects of radiation traversing the human body. In 
particular, radiation has been used in cancer therapy for a long time, and recently this has 
been extended to include heavy ion particle beams. The advent of this particular type of 
therapy has led to the need for increased capabilities in the transport codes used to 
simulate the detailed nature of the treatment doses to the various tissues that are 
encountered. This capability is also of interest to NASA because of the nature of the 
radiation environment in space. [1] While in space, the crew members’ bodies are 
continually being traversed by virtually all forms of radiation. In assessing the risk that 
this exposure causes, heavy ions are of primary importance. These arise both from the 
primary external space radiation itself, as well as fragments that result from interactions 
during the traversal of that radiation through any intervening material including 
intervening body tissue itself. Thus the capability to characterize the details of the 
radiation field accurately within a human body subjected to such external “beams” is of 
critical importance. 

In order to provide this capability by making use of the widest practical 
application of the known physics, the FLUKA Monte Carlo code has been extended over 
the last several years to include the ability to simulate heavy-ion interactions more 
completely. [2] Currently, FLUKA is available with internal event generators that are 
capable of simulating inelastic nuclear reactions down to an incident lab kinetic energy of 
-100 MeV/A.[3] Work is in progress to extend the internal event generator capability all 


the way down to the reaction thresholds, as well as to update and improve the existing 
capabilities on a continuing basis as new data become available. 

One common hurdle in the use of all Monte Carlo radiation transport codes is the 
difficulty of providing the detailed geometry information about the system to be modeled. 
This is particularly problematic when the object is as complex as the human body. This 
situation is multiplied in difficulty by the requirement not only for an accurate positional 
geometric description, but also the absolute need to attach to that positional information 
the details of the composition and density of the actual material. In cases where the 
human body is the object of the simulation, the technique of using information from CT- 
scans has been common, especially where the nature of the incident beam is restricted to 
electromagnetic radiation and electrons. This is because the primary information from a 
CT-scan is a measure of the local distribution of electron density in the object scanned. 
Raw CT-scan files are generally represented as scan layers of fixed thicknesses with data 
from individual square cells within each sequential layer being given in raster pixel 
format. The three-dimensional volume element made up of each of these pixels in a 
given layer is referred to as a “voxel.” 

In order to be able to accept these raw CT-scan inputs, FLUKA has been modified 
to allow for the direct embedding of a region composed of these internal voxels within 
any normal FLUKA geometric description. So for example, one could embed a CT-scan- 
based voxel human phantom region within a spacecraft that was otherwise described 
using the normal FLUKA capabilities. Similarly, a CT-scan-based voxel phantom can be 
embedded in the external laboratory environment geometry description for simulation of 
accelerator-based exposures. 



Unfortunately, by itself, a raw CT-scan is not as useful in providing a basis for the 
transport geometry when hadronic beams are applied due to the general lack of direct 
compositional information. This can be remedied if there is some external process 
whereby the composition can be overlaid. In some cases, like the soft tissue in the human 
body, the use of a common global generic composition is a very good approximation. 
However, when dealing with the whole body or large fractions of it this technique 
produces incorrect results because some organs, particularly bone and the lungs have 
both notably different compositions and porous structures. 

ASSIGNING COMPOSITION 

In order to address these issues, we have explored the development of an 
algorithm to distinguish bone and lung tissue from the other soft tissues within the body. 
The additional small problem of the existence of gas pockets, for example in the stomach 
and bowel, are easily included due to the very dramatic difference in density. 

Our intention is to provide a compositional tag for each voxel that is either a 
generic “soft tissue,” that is global for all soft tissue, or one of the special tissue types 
including hard bone, trabecular bone and porous lung. It should be emphasized that the 
actual measured densities will be used scaled to the generic composition assigned based 
on the x-ray cross sections, which are typically proportional to the electron density. For 
example, the density of each soft tissue voxel can still vary, even though the presumed 
elemental compositions remains the same. 



The challenges for bone and lung tissue are actually of two distinct kinds. First, 
there is the initial problem of simply distinguishing this tissue from the surrounding soft 
tissue. This is generally a pattern recognition problem which yields to well-know 
techniques that take advantage of density variations. Lung tissue is relatively easy to 
distinguish, but bone is more problematic. The thin “hard bone” shell that surrounds 
most human bone tissue is filled with an inner trabecular bone that is a combination of 
hard-bone filaments and a variety of different kinds of marrow. This composite 
trabecular bone tissue is typically unresolved on many large-scale CT-Scans, and has an 
average density that can be similar to the surrounding soft tissue. However, the general 
property that it is internal to the hard-bone shell typically provides a handle on 
distinguishing it from the external soft-tissue. We have developed several software tools 
to facilitate the development of these algorithms. 

In general the work reported here has been done using two separate sources for 
CT-scans. The first is a whole -body image of an adult male [2] that has been previously 
analyzed by hand to assign voxel by voxel the organ membership along with the 
corresponding material composition. The second image is that of trabecular bone and 
represents a 1.2 cm cubic volume entirely with the region of trabecular bone with a 40 
pm voxel size. In this case, all voxels are identified either as entirely hard-bone or 
entirely interstitial marrow. [4] 

Figure 1 shows the displays employed to develop the appropriate algorithms. 
Figure 1(a) contains the whole body image where the voxel densities of the slices shown 
have been displayed as a corresponding linear gray-scale. Figure 1(b) shows the same 
image processed by a simple density-based algorithm to distinguish the hard bone from 



soft tissue. Note that there is a histogram tool towards the top which displays the actual 
density profile for the vertical cross-hair line and a fit to that profile by the algorithm. 
The results are displayed for comparison in the image itself. One can see that the hard 
bone is easily distinguished from the soft tissue in a general sense. The major difficulties 
are with the boundary voxels and distinguishing the internal trabecular bone from general 
soft tissue. 

The boundary between hard bone and soft tissue can be problematic because the 
actual physical boundary between the outer hard-bone surface and the surrounding soft 
tissue generally occurs in the midst of a voxel rather than close to or along a voxel side. 
The consequence is a density for the voxel that is intermediate between that of the soft 
tissue and the hard bone. One possibility is to treat this as a special material that is a 
weighted homogeneous mixture of bone and soft tissue with the weighting being a linear 
fit based on the intermediate density value itself. This has the seeming benefit of 
including the effects of the nearby bone or soft tissue, whichever is the minority, in 
contrast with the our current algorithm which simply uses a mean value as a threshold 
either to declare the entire voxel as being bone or soft tissue. However, the actual result 
is to guarantee that almost all tracks traversing the voxel will have to deal with an 
incorrect mix, as only a small number will actually see the same medium as the average. 
On the other hand, our algorithm has the benefit that at least some of the tracks will be 
seeing the correct density. This is the simplest reasonable solution. A more complex 
solution that may be employed is similar to the one that has been developed to deal with 
the trabecular bone regions as discussed in the next section. 


Figure 1(b) also shows that the lungs are easily distinguished from the 
surrounding tissue, and even the central region where the major pulmonary arteries and 
veins enter the lungs, the distinctions are fairly clear. The lung boundaries present the 
same challenge as does the bone boundary, and are amenable to the same solutions. 
These simple algorithms are subject to improvement, but given the level of success that 
has already been achieved by our methods, considerable progress has been made in 
getting the material properties into the Monte Carlo geometry inputs. 


TRABECULAR BONE POROSITY 


The porosity of trabecular bone and the lungs provides its own challenge. In this 
case, the problem stems from the fact that the details of porosity occur at a scale that is 
not practical to evaluate in large volume whole-body CT-scans. As such if the voxel is 
treated as possessing some mean density with a homogeneous composition, then all 
tracks traversing the voxel with the same pathlength in the voxel will see the identical 
medium. In reality, the pseudo-random porous nature of the material will actually give 
rise to a distribution of pathlengths. As we shall show, this distribution can be very 
directionally dependent in the case of trabecular bone. 

Employing the CT-scan image of trabecular bone with 40 pm resolution 
mentioned in the previous section, several FLUKA runs were performed using rays to 
determine the actual pathlengths within the hard bone component of the trabecular bone. 
In FLUKA, rays are imaginary particles that travel in straight lines and can be used to 
account for the details of the material traversed along their trajectories. This particular 



example of trabecular bone possessed about 10% of its voxels as hard-bone with the 
remainder being treated as vacuum for our purposes here. Hard bone is taken to have a 
density of 1.920 gm/cm 2 and the pathlength distributions are all given in terms of the 
gm/cm 2 of hard bone traversed per cm of total pathlength traversed. The data are limited 
to rays with a minimum 0.5 cm pathlength in the trabecular bone volume. 

Figure 2(a) shows the pathlength distribution for -100,000 rays which are 
incident isotropically on the CT-scan volume. It is seen that there are essentially no rays 
that manage to avoid all hard-bone voxels, and the mean is about 0.1 gm/cm /cm. Figure 
2(b) is an expanded view of the long pathlength tail of the same distribution. 

Figure 3(a) shows the result for another 100,000 rays, but for which the direction 
is strictly vertical, along the axis of the bone. Slightly more than 50% of these rays pass 
through the full 1.2 cm of trabecular bone without striking any hard bone voxels!.. This 
implies that strictly along the axis of the bone, the pathlength variation is a very 
important feature for particle transport. Figure 3(a) shows a similar set of rays, but where 
the angle with respect to the bone axis is distributed uniformly over a 100 mrad cone. 
However, note that even this modest angular deviation from the strict longitudinal axis 
causes a rapid reduction of this anomaly with the result tending towards the isotropic 
distribution with the virtual disappearance of the zero pathlength rays. 

Figure 4 shows the result for a lateral beam (perpendicular to the bone axis) with 
no angular spread. A comparison with Figure 2(a) shows that the lateral pathlength 
distribution is comparable to the general isotropic distribution. 

Thus, it appears as if a strategy can be deployed wherein the pathlength used is 
compensated by a normal distribution fit to the isotropic pathlength distribution, with the 



provision that for particle directions that are very close to axial along the bone, a second 
distribution should be used that is a combination of (a) 50% no hard bone, and (b) 50% of 
twice the standard isotropic pathlength distribution. 

To implement this strategy, rather than vary the actual pathlength, the mechanism 
to be employed within FLUKA will be to provide for a new kind of material that has a 
variable density. Then when a particle enters such a medium the program will select the 
density to be used for that particle randomly from a distribution function. This will allow 
the effect of the porosity of the medium to be properly simulated and the net effect on a 
traversing beam of charged particles will be to spread out the energy loss effects as well 
as the other interactions probabilities. 

CONCLUSIONS 


FLUKA has the capability to employ a voxel geometry, and that feature can be 
used to import information directly from raw CT-scan output files. We have developed a 
simple algorithm to classify tissue into 4 generic categories. These include soft-tissue, 
hard bone, trabecular bone and porous lung. With that categorization we can assign 
generic compositions, and actual densities using the measured CT-scan values to scale 
the density voxel by voxel. Finally, we have explored the structure of trabecular bone in 
some detail and have determined the pathlength distribution of typical rays in hard bone 
within that medium. We have highlighted the substantial asymmetry in the distribution 
as the rays tend to be aligned along the bone’s longitudinal axis. It is planned to 
implement a new type of material within the FLUKA geometry that will allow a variable 



density. Such a variable density material will be able to employ a value that is selected 
from a specified distribution on a track by track basis. A similar approach is anticipated 
for the treatment of porous lung tissue, but the detailed examination of the pathlength 
distributions in lung structure is still being explored. 

ACKNOWLEDGMENTS 


We would like to acknowledge the use of the so-called “Golem” whole-body CT- 
scan data set along with the version that resulted from the efforts of Laura De Biaggi in 
her detailed voxel by voxel assignment of individual organs to that data set. We are also 
grateful to Professor Gemunu Gunaratne of the University of Houston for providing 
access to the trabecular bone CT-scan employed as well as for his efforts in producing 
graphic images of the CT-scan. Finally, we gratefully acknowledge the permission of 
Professor Michael A.K. Liebschner of the Bioengineering Department of Rice University 
to use the trabecular bone CT-scan which Professor Gunaratne had obtained from him 
originally. This work was partially supported by the EC (contract no. FI6R-CT-2003- 
508842, "RISC-RAD"), NASA (Contracts NAG8-1901 and NAG8-1658) and the 
Institute for Space Systems Operations at the University of Houston. 

REFERENCES 

1. V. Andersen F. Ballarini G. Battistoni, M. Campanella, M. Carboni, F. Cerutti, 
A. F.m pl , A. Fasso, A. Ferrari, E. Gadioli, M.V. Garzelli, K. Lee, A. Ottolenghi, 


M. Pelliccioni, L.S. Pinsky, J. Ranft, S. Roesler, P.R. Sala and T.L. Wilson, The 
FLUKA code for space applications: recent developments. Adv. Space Res, in 
press (2004). 

2. F. Ballarini, M. Biaggi L. De Biaggi, A. Ferrari, A. Ottolenghi, A. Panzarasa, 
Ft.G. Paretzke, M. Pelliccioni, P. Sala, D. Scannicchio and M. Zankl, Role of 
shielding in modulating the effects of Solar Particle Events: Monte Carlo 
calculation of physical and “biological” dose in different organs. Adv. Space 
Res., in press (2004). 

3. Fasso, A. Ferrari, S. Roesler, P.R. Sala, F. Ballarini, A. Ottolenghi, G. Battistoni, 
F. Cerutti, E. Gadloli, M.V. Garzelli, A. Empl, J. Ranft The physics models of 
FLUKA: status and recent developments CHEP-2003-MOMT005, Jun 2003. 
10pp. Talk given at 2003 Conference for Computing in High-Energy and Nuclear 
Physics (CHEP 03), La Jolla, California, 24-28 Mar 2003. Published in eConf 
C0303241 :MOMT005,2003 e-Print Archive: hep-ph/0306267 (2003). 


4. Gunaratne, G. and Leibschner, M., Private Communication (2004). 



FIGURE CAPTIONS 


Figure 1. Screen Shots are presented of the software tool developed to assess the 
performance of the tissue assignment algorithm. Figure 1(a) shows a slice from the 
“Golem” CT-scan where the raw density values are displayed on a linear brightness 
scale as pixels. Figure 1(b) shows the result of filtering this image with the simple 
threshold algorithm. Note the histogram tool that displays the raw data across the 
slice and the values from the algorithm superimposed. The Z-axis is the vertical body 
axis and the X-axis is lateral from side to side, with the Y-axis being front to back. 
Positive Y is in the forward direction from the body’s face. 

Figure 2. This displays the pathlength distributions in trabecular bone for an isotropic 

irradiation. All pathlengths included in the plot have a minimum 0.5 cm absolute 

length, and the results are plotted as gm/cm 2 of hard bone per actual cm of pathlength. 

Figure 2(a) is an enlargement of the long pathlength tail. The distribution fits a 

normal distribution quite well. 

* 

Figure 3. This displays the pathlength distributions similar to those of Figure 1, but 
for tracks incident along the longitudinal axis of the bone (the Z-axis) in our 
coordinate system. Figure 2(a), which is a plot for rays strictly parallel to the Z-axis, 
shows the surprising behavior of possessing a considerable peak at zero pathlength. 
However, in contrast, Figure 2(b), which has its rays diverging from the Z-axis with 
this divergence spread uniformly over angle up to 100 mRadians, and the large zero 


pathlength peak is entirely missing. 


Figure 4. This plot shows the pathlength distribution that results from an exposure to 

* 

rays that are exactly parallel to one of the axes that is perpendicular (the reference to 
X-axis here is intended to be a generic axis perpendicular to the longitudinal axis of 
the bone). Like Figure 3(a), there is no divergence in the beam, but unlike that figure, 
there is no corresponding zero pathlength peak. 
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Figure 2(a) 


Long Pathlength Distribution Tail for Isotropic Rays in Trabecular Bone 







# Per 0.01 gm/cm A 2/cm Bin # Per 0.020 gm/cm A 2/cm Bin 


Pathlength Distribution for Vertical Rays in gm/cm A 2/cm 


50000 - 


40000 - 


30000 — 


20000 - 


10000—1 


Statistics 

Entries 

100000 1 

Mean 

-'BEE 

RMS 

“11 


0.05 0.1 0.15 0.2 0.25 0.3 0.35 

Pathlength in gm/cm A 2/cm 


0.4 0.45 0.5 


Figure 3(a) 


Pathlength Distribution for Z-Axis Beam With 0.1 mRad Divergence 



Figure 3(b) 
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